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Abstract 

The Large Area Telescope (LAT), the main instrument of the Fermi Gamma-Ray Space 
Telescope, detects high energy gamma rays with energies from 20 MeV to more than 300 GeV. 
The two main scientific objectives, the study of the Milky Way diffuse background and the de- 
tection of point sources, are comphcated by the lack of photons. That is why we need a powerful 
Poisson noise removal method on the sphere which is efficient on low count Poisson data. This 
paper presents a new multiscale decomposition on the sphere for data with Poisson noise, called 
Multi-Scale Variance Stabilizing Transform on the Sphere (MS-VSTS). This method is based 
on a Variance Stabilizing Transform (VST), a transform which aims to stabilize a Poisson data 
set such that each stabilized sample has a quasi constant variance. In addition, for the VST 
used in the method, the transformed data are asymptotically Gaussian. MS-VSTS consists of 
decomposing the data into a sparse multi-scale dictionary like wavelets or curvelets, and then 
applying a VST on the coefficients in order to get almost Gaussian stabilized coefficients. In 
this work, we use the Isotropic Undecimated Wavelet Transform (lUWT) and the Curvelet 
Transform as spherical multi-scale transforms. Then, binary hypothesis testing is carried out 
to detect significant coefficients, and the denoised image is reconstructed with an iterative al- 
gorithm based on Hybrid Steepest Descent (HSD). To detect point sources, we have to extract 
the Galactic diffuse background: an extension of the method to background separation is then 
proposed. In contrary, to study the Milky Way diffuse background, we remove point sources 
with a binary mask. The gaps have to be interpolated: an extension to inpainting is then pro- 
posed. The method, applied on simulated Fermi LAT data, proves to be adaptive, fast and easy 
to implement. 
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1. Introduction 

The Fermi Gamma-ray Space Telescope, 
which was launched by NASA in June 2008, 
is a powerful space observatory which studies 
the h igh-energy gamma-ray sky (jAtwood et aLl 
120091 ^. Fermi's main instrument, the Large Area 
Telescope (LAT), detects photons in an energy 
range between 20 MeV to greater than 300 GeV. 
The LAT is much more sensitive than its prede- 
cessor, the EGRET telescope on the Compton 
Gamma Ray Observatory, and is expected to 
find several thousand gamma ray sources, which 
is an order of magnitude mor e than its predeces- 
sor EGRET (jHartman et al.lflQQa) . 

Even with its ~ Im^ effective area, the num- 
ber of photons detected by the LAT outside the 
Galactic plane and away from intense sources 
is expected to be low. Consequently, the spher- 



ical photon count images obtained by Fermi 
are degraded by the fluctuations on the number 
of detected photons. The basic photon-imaging 
model assumes that the number of detected pho- 
tons at each pixel location is Poisson distributed. 
More specifically, the image is considered as a 
realization of an inhomogeneous Poisson pro- 
cess. This quantum noise makes the source de- 
tection more difficult, consequently it is better 
to have an efficient denoising method for spher- 
ical Poisson data. 

Several techniques have been proposed in 
the literature to estimate Poisson intensity in 
2D. A major class of methods adopt a mul- 
tiscale bayesian framework specifically tailored 
for Poisson data (Novralc&JKolaczyk 2000), in- 
depeii dentl y initiated by iTirnmerman fc Nowakl 
(I1S99D and lKolaczvd (119991) . I Lefkimmiaits et aLl 
(2009) proposed an improved bayesian frame- 
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work for analyzing Poisson processes, based on a 
multiscale representation of the Poisson process 
in which the ratios of the underlying Poisson 
intensities in adjacent scales are modeled as 
mixtures of conjugate parametric distributions. 
Another approach includes preprocessing the 
count data by a va riance stabi l izing transform 
(VST ) such as the lAnscombd (|l948f ) and the 
iFisd (|l955) transforms, applied respectively in 
the spat ial iD onoho 1993) or in the wavelet 
domain (jFrvzlewicz fc Nasonll2004[ ). The trans- 
form reforms the data so that the noise ap- 
proximately becomes Gaussian with a constant 
variance. Standard techniques for independant 
identically distrib uted Gaussian noise are then 
used for denoising. I Zhang et al.l (j2008l) proposed 
a powerful method called Multi-Scale Variance 
Stabilizing Tranform (MS- VST). It consists in 
combining a VST with a multiscale trans- 
form (wavelets, ridgelets or curvelets), yielding 
asymptotically normally distributed coefficients 
with known variances. The choice of the multi- 
scale method depends on the morphology of the 
data. Wavelets represent more efficiently regular 
structures and isotropic singularities, whereas 
ridgelets are designed to represent global lines 
in an image, and curvelets represent efficiently 
curvilinear contours. Significant coefficients are 
then detected with binary hypothesis testing, 
and the final estima te is reconstructed with an 
iterative scheme. In IStarck et al.l (120091 ) . it was 
shown that sources can be detected in 3D LAT 
data (2D+time or 2D+energy) using a specific 
3D extension of the MS- VST. 

There is, to our knowledge, no method for 
Poisson intensity estimation on spherical data. 
It is possible to decompose the spherical data 
into several 2D projections, denoise each pro- 
jection and reconstitute the deiioised spherical 
data, but the projection induces some caveats 
like visual artifacts on the borders or deforma- 
tion of the sources. 

In the scope of the Fermi mission, we have 
two main scientific objectives: 

— Detection of point sources to build the cata- 
log of gamma ray sources, 

— Study of the Milky Way diffuse background, 
which is due to interaction between cosmic 
rays and interstellar gas and radiation. 

The first objective implies the extraction of 
the galactic diffuse background. Consequently, 
we want a method to suppress Poisson noise 
while extracting a model of the diffuse back- 
ground. The second objective implies the sup- 
pression of the point sources: we want to apply 
a binary mask on the data (equal to on point 
sources, and to 1 everywhere else) and to denoise 
the data while interpolating the missing part. 
Both objectives are linked: a better knowledge 
of the Milky Way diffuse background enables us 



to improve our background model, which leads 
to a better source detection, while the detected 
sources are masked to study the diffuse back- 
ground. 

The aim of this paper is to introduce a 
Poisson denoising method on the sphere called 
Multi-Scale Variance Stabilizing Transform on 
the Sphere (MS-VSTS) in order to denoise 
the Fermi photon count maps. This method 
is based on the MS- VST (jZhang et al.l I2008D 
and on recent on mu l ti-sca. l e transforms on the 
sphere (Starck et all l2006t lAbrial et all 120071: 
Abrjal et al, 20081) . Section 2 recalls the multi- 
scale transforms on the sphere which are used 
in this paper, and Gaussian denoising methods 
based on sparse representations. Section 3 intro- 
duces the MS-VSTS. Section 4 applies the MS- 
VSTS to spherical data restoration. Section 5 
applies the MS-VSTS to inpainting. Section 6 
applies the MS-VSTS to background extrac- 
tion. Conclusions are drawn in Section 7. In 
this paper, all experiments were performed on 
HEA LPix maps with nside = 128 (jGorski et al.l 
|2005( ). which corresponds to a good pixelisation 
choice for the GL AST/FERMI resolution. The 
performance of the method is not dependent 
on the nside parameter. For a given data set, 
if nside is small, it just means that we don't 
want to investigate the finest scales. If nside is 
large, the number of counts per pixel will be 
very small, and we may not have enough statis- 
tics to get any information at the finest reso- 
lution levels. But it will not have any bad ef- 
fect on the solution. Indeed, the finest scales 
will be smoothed, since our algorithm will not 
detect any significant wavelet coefficients in the 
finest scales. Hence, starting with a fine pixeli- 
sation (i.e. large nside), our method will pro- 
vide a kind of automatic binning, by threshold- 
ing wavelets coefficients at scales and at spatial 
positions where the number of counts is not suf- 
ficient. 



2. Multi-scale analysis on the sphere 

New multi-scale transforms o n the sphere 
were developed by IStarck et ahl (|2006[ ). These 
transforms can be inverted and are easy to com- 
pute with the HEALPix pixellisation, and were 
used for denoising, deconvolution, morphologi- 
cal co mponent ana lysis and impainting applica- 
tions (jAbrial et al.ll2007f ). In this paper, we use 
the Isotropic Undecimated Wavelet Transform 
(lUWT) and the Curvelet Transform. 
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2.1. Multi-scale Transforms on the sphere 

2.1.1. Isotropic Undecimated Wavelet Transform 
on the sphere 

The Isotropic Undecimated Wavelet 
Transform on the sphere (lUWT) is a wavelet 
transform on the sphere based on the spherical 
harmonics transform and with a very simple 
reconstruction algorithm. At scale j, we denote 
aj{6,ip) the scale coefRcients, and dj(6,ip) 
the wavelet coefficients, with 9 denoting the 
longitude and ip the latitude. Given a scale 
coefficient aj, the smooth coefficient flj+i is 
obtained by a convolution with a low pass filter 
hj : flj+i — ttj * hj. The wavelet coefficients 
are defined by the difference between two 
consecutive resolutions : dj+i = Qj — Oj+i. A 
straightforward reconstruction is then given by: 

J 

ao{e,^)^aj{e,^)+Y,djiO,^) (1) 
i=i 

Since this transform is redundant, the procedure 
for reconstructing an image from its coefficients 
is not unique and this can be profitably used to 
impose additional constraints on the synthesis 
functions (e.g. smoothness, positivity). A recon- 
struction algorithm based on a variety o f filter 
banks is described in IStarck et al.l (|2006( ) . 

2.1.2. Curvelet Transform on the sphere 

The curvelet transform enables the direc- 
tional analysis of an image in different scales. 
The data undergo an Isotropic Undecimated 
Wavelet Transform on the sphere. Each scale 
j is then decomposed into smoothly overlap- 
ping blocks of side-length Bj in such a way 
that the overlap between two vertically adjacent 
blocks is a rectangular array of size Bj x Bj/2, 
using the HEALPi x pixellisation. Finally, the 
ridgelet transform (jCandes &: Donohol |1999[ ) is 
applied on each individual block. The method 
is best for the detection of anisotropic struc- 
tures and smooth curves and edges of dif- 
ferent lengths. Mo re details can be found in 
IStarck et al.l (|20nfiD . 

2.2. Application to gaussian denoising on the 
sphere 

Multiscale transforms on the sphere have 
been used successfully for Gaussian denoising 
via non-linear filtering or thresholding methods. 
Hard thresholding, for instance, consists of set- 
ting all insignificant coefRcients (i.e. coefficients 
with an absolute value below a given thresh- 
old) to zero. In practice, we need to estimate 
the noise standard deviation aj in each band j 
and a coefficient Wj is significant if \wj\ > nuj, 



where k is a parameter typically chosen between 
3 and 5. Denoting Y the noisy data and HT\ 
the thresholding operator, the filtered data X 
are obtained by: 

X = *i?rA(*^Y), (2) 

where 4>-^ is the multiscale transform (lUWT or 
curvelet) and $ is the multiscale reconstruction. 
A is a vector which has the size of the number 
of bands in the used multiscale transform. The 
thresholding operation thresholds all coefficients 
in band j with the threshold \j = naj . 

3. Multi-Scale Variance Stabilizing 
Transform on the Sphere (MS-VSTS) 

3.1. Principle of VST 

3.1.1. VST of a Poisson process 

Given Poisson data Y := {Yi)i, each sam- 
ple Yi ~ 'P{Xi) has a variance Var[yi] = A^. 
Thus, the variance of Y is signal-dependant. 
The aim of a VST T is to stabilize the data such 
that each coefficient of T(Y) has an (asymptot- 
ically) constant variance, say 1, irrespective of 
the value of A^. In addition, for the VST used in 
this study, T(Y) is asymptotically normally dis- 
tributed. Thus, the VST-transformed data are 
asympt otically stationary and gaussian. 

The lAnscomb^ (|1948( ) transform is a widely 
used VST which has a simple square-root form 

T{Y) := 2v/r + 3/8. (3) 

We can show that T(Y) is asymptotically nor- 
mal as the intensity increases. 

T{Y) - 2VA >Af{0, 1) (4) 

A —?' +00 

It can be shown that the Anscombe VST re- 
quires a high underlying intensity to well stabi- 
lize t he data (typically for A ^ 10) l|Zhang et al.l 
I2008D . 

3.1.2. VST of a filtered Poisson process 

Let Zj :— J2i be the filtered process 

obtained by convolving {Yi)i with a discrete fil- 
ter h. We will use Z to denote any of the Zj's. 
Let us define := ^'^^ ^ — 1? 2, • • • . 

In addition, we adopt a local homogeneity as- 
sumption stating that Xj-i — A for all i within 
the support of h. 

We define the square-root transform T as fol- 
lows: 

T{Z):=b-aigniZ + c)\Z + c\^/^, (5) 

where & is a normalizing factor. Lemma [T] proves 
that T is a VST for a filtered Poisson pro- 
cess (with a nonzero-mean filter) in that T(Y) 
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is asymptotically normally distributed with a 
stabilized variance as A becomes large (see 
IZhang et all (|2008l ) for a proof). 

Lemma 1 (Square root as VST) If ti ^ 0, 

II^IU, II^IU < oo, then we have : 



sign(Z + c)y/\Z + c\ - sign(Ti) Vln|A 



V 



X — > +00 



-MO 



T2 

4|ri| 



(6) 



3.2. MS-VSTS 

The MS-VSTS consists in combining the 
square-root VST with a multi-scale transform. 



3.2.1. MS-VSTS + lUWT 

This section describes the MS-VSTS + 
lUWT, which is a combination of a square-root 
VST with the lUWT. The recursive scheme is: 



lUWT 



2j = * Oj-l 



MS-VSTS r aj = hj^i * a,_i 
' + lUWT 1 d, = rj-i(aj-i) - T,{aj) 



(7) 



In ([7]), the filtering on flj-i can be rewritten 
as a filtering on ao '= Y, i.e., aj ~ /i^-'-' * oq, 
where /i'^' = ^j-i * • • • * ^1 * /io for j ^ 1 and 
= S, where S is the Dirac pulse (J = 1 on 
a single pixel and everywhere else). Tj is the 
VST operator at scale j: 



T,{a,) = 6(^")sign(a, + c^^'^)^ \a, + cU)\. (8) 



Let us define t. 



In IZhang et al.l ()2008l ). it has ben shown 
that, to have an optimal convergence rate for 
the VST, the constant c^-'^ associated to h^^^ 
should be set to: 



7t, 



U) 



Ao) 



U) 



2t. 



(i) 



(9) 



The MS-VSTS+IUWT procedure is directly in- 
vertible as we have: 



(0,^). (10) 



Setting b'^^^ := sgn(T}-'^ ) / \J \t[^^ | , if A is constant 
within the support of the filter. h'^j\ then we 



have (jZhang et al.ll2008[ ): 
V ( 



dj{e,if)- 



►AA 0, 



A ^ -l-oo 



4t, 



0-1)" 



(11) 



4r| 



2rl 



where (., .) denotes inner product. 

It means that the detail coefficients issued 
from locally homogeneous parts of the signal 
follow asymptotically a central normal distri- 
bution with an intensity-independant variance 
which relies solely on the filter h and the current 
scale for a given filter h. Consequently, the sta- 
bilized variances and the constants b^^\S-^\T^'^ 
can all be pre-computed. Let us define a^-^ the 
stabilized variance at scale j for a locally homo- 
geneous part of the signal: 



.0) 



(/i(j-i),/i«) 



4Tj' 



2rl 



(12) 



To compute the ctq), b^^\c^^\Tj,''\ we only 

have to know the filters h^^\ We compute these 
filters thanks to the formula aj ~ h^^'^ * op, by 
applying the lUWT to a Dirac pulse oq = 6. 
Then, the /i^^^ are the scaling coefficients of the 
lUWT. The (T(j) have been precomputed for a 
6-scaled lUWT (Tabled]). 

We have simulated Poisson images of differ- 
ent constant intensities A, computed the lUWT 
with MS-VSTS on each image and observed 
the variation of the normalized value of (T(j) 

((0'(j))simulated/ (o-(j) )theoretical) aS a fuUCtiou of A 

for each scale j (Fig.[T]). We see that the wavelet 
coefficients are stabilized when A > 0.1 except 
for the first wavelet scale, which is mostly con- 
stituted of noise. On Fig. [51 we compare the 
result of MS-VSTS with Anscombe + wavelet 
shrinkage, on sources of varying intensities. We 
see that MS-VSTS works well on sources of very 
low intensities, whereas Anscombe doesn't work 
when the intensity is too low. 

3.2.2. MS-VSTS + Curvelets 

As the first step of the algorithm is an 
lUWT, we can stabilize each resolution level as 
in Equation ([7]) . We then apply the local ridgelet 
transform on each stabilized wavelet band. 

It is not as straightforward as with the 
lUWT to derive the asymptotic noise variance 
in the stabilized curvelet domain. In our experi- 
ments, we derived them using simulated Poisson 
data of stationary intensity level A. After hav- 
ing checked that the standard deviation in the 
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Figure 1. Normalized value ((o'(j))siinuiated/(o'(j))thcoreticai) of the stabilized variances at each scale 
j as a function of A. 



curvelet bands becomes stabilized as the inten- 
sity level increases (which means that the stabi- 
lization is working properly), we stored the stan- 
dard deviation aj^i for each wavelet scale j and 
each ridgelet band I (Table [5]) . 

4. Poisson Denoising 

4.1. MS-VST + lUWT 

Under the hypothesis of homogeneous 
Poisson intensity, the stabilized wavelet coeffi- 
cients dj behave like centered Gaussian variables 
of standard deviation a(j) . We can detect signif- 
icant coefficients with binary hypothesis testing 
as in Gaussian denoising. 

Under the null hypothesis Ho of homoge- 
neous Poisson intensity, the distribution of the 
stabilized wavelet coefficient dj [k] at scale j and 
location index k can be written as: 



p{d,[k]) = 



1 



2Tra. 



■ exp{-dj[k]^ /2a'j 



(13) 



The rejection of the hypothesis Hq depends 
on the double-sided p-value: 



Pj[k] 



1 



+ 00 



lrf,[fe]l 



exp(— a; /2a^)dx. (14) 



Consequently, to accept or reject Hq, we 
compare each with a critical threshold 

KtTj, K — 3, 4 or 5 corresponding respectively 
to significance levels. This amounts to deciding 
that: 

— if ^ naj, dj[k] is significant. 

— if < KCTj, dj[k] is not significant. 

Then we have to invert the MS-VSTS scheme 
to reconstruct the estimate. However, although 
the direct inversion is possible (Eq. (fTU|)). it 
can not guarantee a positive intensity estimate, 
while the Poisson intensity is always nonnega- 
tive. A positivity projection can be applied, but 
important structures could be lost in the esti- 
mate. To tackle this problem, we reformulate the 
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reconstruction as a convex optimisation problem 
and solve it iteratively with an algorith m based 
on H ybrid Steepest Descent (HSD) f|Yamadal 
I2001h . 

We define the multiresolution support A^, 
which is determined by the set of detected sig- 
nificant coefficients after hypothesis testing: 

Al := {{j, fc)|if dj[k] is declared significant}. 

(15) 

We formulate the reconstruction problem as 
a convex constrained minimization problem: 



Arg min II* X||i,s.t. 



V(j, fc)GX,(*^X),[fc]-(*^Y),[fc] 



Algorithm 1 MS-VSTS + lUWT Denoising 
Require: data ao 

tions iVmax, threshold k 

Detection 



Y, number of itera- 



(16) 



5 
6 
7 

8 

9 

10: 
11: 



for j = 1 to J do 
Compute flj and d 
Hard threshold 
and update Ai. 

end for 

Estimation 



J using (O. 
dj [k] I with threshold naj 



Initialize X(o) = 0, Aq = 1. 
for n = to A'max — 1 do 

X = P+[X(") + *Pm*^(Y - X("))]. 

X("+i) = *STa,.[*'^X]. 

"+1 ~ — N^„^-l — • 

end for 

Get the estimate A = X^^"--). 



where # denotes the lUWT synthesis operator. 

This problem is solved with the following it- 
erative scheme: the image is initialised by X^^^ = 
0, and the iteration scheme is, for n = to 
A^max - 1: 



X = F+[X(") +*Fm*^(Y-X("))] (17) 
X(«+i) = *STaJ*^X] (18) 



At each iteration n, we compute the MS- 
VSTS of X("). We denote d^f^k] the stabilised 

coefficients of X*^"-' . We make a hard threshold- 
ing on {dj[k]-d'-^\k]) with the same thresholds 
as in the detection step. Significant coefficients 
are added to the multiresolution support A4. 



where P+ denotes the projection on the posi- 
tive orthant, Pm denotes the projection on the 
multiresolution support Ai: 



PMdj[k] 



dj[k] if {j,k)&M, 




otherwise 



(19) 



and STa„ the soft-thresholding with threshold 



_ r sign(d)(|d| -A„) ifM|^A„, 







otherwise 



(20) 

We chose a decreasing threshold A„ — 

N„ 



N„- 



The final estimate of the Poisson intensity 
is: A = X'^™"''^. Algorithm [T] summarizes the 
main steps of the MS-VSTS + lUWT denoising 
algorithm. 

4.2. Multi-resolution support adaptation 

When two sources are too close, the less in- 
tense source may not be detected because of 
the negative wavelet coefficients of the bright- 
est source. To avoid such a drawback, we may 
update the multi-resolution support at each it- 
eration. The idea is to withdraw the detected 
sources and to make a detection on the remain- 
ing residual, so as to detect the sources which 
may have been missed at the first detection. 



Algorithm 2 MS-VSTS + lUWT Denoising 
Multiresolution Support Adaptation 



Require: data gq :— Y, 
tions A'^inax, threshold k 
Detection 



number of itera- 



5 
6 
7 
8 
9 

10: 

11 
12 
13 



for j = 1 to J do 

Compute aj and dj using ([7|. 

Hard threshold |rfj[fc]| with threshold naj 

and update A4. 
end for 
Estimation 



Initialize X'"' 0, Ao 1. 
for n = to A'^inax — 1 do 

X = P+[X(") + *Pa^*^(Y - X("))]. 

X("+i) = $STa„[*^X]. 

Compute the MS-VSTS on X(") to get the 



stabilised coefFcients d 



Hard threshold |c?j[/c] — (ij"-'[fc]| and update 
M. 

A,i+i = 
end for 

Get the estimate A = xt"""-). 



-("+!) 



N„ 



The main steps of the algorithm are sum- 
marized in Algorithm [51 In practice, we use 
Algorithm [2] instead of Algorithm [1] in our ex- 
periments. 
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4.3. MS-VST + Curvelets 

Insignificant coefiicients are zeroed by us- 
ing the same hypothesis testing framework as 
in the wavelet scale. At each wavelet scale j 
and ridgelet band fc, we make a hard threshold- 
ing on curvelet coefficients with threshold Kaj,k, 
K — 3, 4 or 5. Finally, a direct reconstruction 
can be performed by first inverting the local 
ridgelet transforms and then inverting the MS- 
VST + lUWT (Equation Uni)). An iterative re- 
construction may also be performed. 

Algorithm [3] summarizes the main steps of 
the MS-VSTS + Curvelets denoising algorithm. 



Algorithm 3 MS-VSTS + Curvelets Denoising 
1: Apply the MS-VST + lUWT with J scales 
to get the stabilized wavelet subbands dj. 

2: Set Bi = Sniin- 

3: for j = 1 to J do 

4: Partition the subband dj with blocks 
of side-length Bj and apply the digital 
ridgelet transform to each block to obtain 
the stabilized curvelets coefficients. 

5: if j modulo 2 = 1 then 

6: Bj+i = 2Bj 

7: else 

8: Sj+i — Bj 

9: end if 

10: HTs on the stabilized curvelet coefficients. 
11: end for 

12: Invert the ridgelet transform in each block 
before inverting the MS-VST + lUWT. 



4.4. Experiments 

The method was tested on simulated Fermi 
data. The simulated data are the sum of a 
Milky Way diffuse background model and 1000 
gamma ray point sources. We based our Galactic 
diffuse emission model intensity on the model 
gll_iem_v02 o btained at t he Fermi Science 
Support Center (jMversI [2009h . This model re- 
sults from a fit of the LAT photons with vari- 
ous gas templates as well as inverse Compton 
in several energy bands. We used a realistic 
point-spread function for the sources, based 
on Monte Carlo simulations of the LAT and 
accelerator tests, that scale approximately as 
0.8{E/lGeV)-°-^ degrees. The position of the 
205 brightest source s were taken from the Fermi 
3-month source list (|Abdo et al.ll2009[) . The po- 
sition of the 795 remaining sourc es follow the 
LAT 1-year Point Source Catalog (|Mversll201(3l) 
sources distribution: each simulated source was 
randomly sorted in a box of Al=5° and Ab=l° 
around a LAT 1-year catalog source. We simu- 
lated each source assuming a power-law depen- 



dence with its spectral index given by the 3- 
month source list and the first year catalog. We 
used an exposure of 3.10^°s.crn^ corresponding 
approximatively to one year of Fermi all-sky sur- 
vey around 1 GeV. The simulated counts map 
shown here correspond to photons energy from 
150 MeV to 20 GeV. 

Fig. E] compares the result of denoising with 
MS-VST + lUWT (Algorithm [J), MS-VST + 
curvelets (Algorithm [3]) and Anscombe VST + 
wavelet shrinkage on a simulated Fermi map. 
Fig. m shows one HEALPix face of the results. 
As expected from theory, the Anscombe method 
produces poor results to denoise Fermi data, 
because the underlyning intensity is too weak. 
Both wavelet and curvelet denoising on the 
sphere perform much better. For this applica- 
tion, wavelets are slightly better than curvelets 

[SNRwavelets = 65 .8dB , SN Rcurvelets = 

37.3dB, SNR{dB) = 201og( 

^signal j ^ noise 

)). As 

this image contains many point sources, thisre- 
sult is expected. Indeed wavelet are better than 
curvelets to represent isotropic objects. 

5. Milky Way diffuse background study: 
denoising and inpainting 

In order to extract from the Fermi photon 
maps the galactic diffuse emission, we want to 
remove the point sources from the Fermi im- 
age. As our HSD algorithm is very close to the 
MCA algorithm (jStarck et al.ll2004l ). an idea is 
to mask the most intense sources and to mod- 
ify our algorithm in order to interpolate through 
the gaps exactly as in the MCA-Inpainting al- 
gorithm (jAbrial et al.l [20071. This modified al- 
gorithm can be called MS-VSTS-Inpainting al- 
gorithm. 

The problem can be reformulated as a convex 
constrained minimization problem: 

Argrnin ||$"^X||i, s.t. 

r x^o, 

\ V(j,fc) e X, (*^nx),[A;] = (*^Y),[fc], 

(21) 

where 11 is a binary mask (1 on valid data and 
on invalid data). 

The iterative scheme can be adapted to cope 
with a binary mask, which gives: 

X = P+[X(") + *Pm *^n(Y - X("))], (22) 
X("+i) ^ #STa„[*X]. (23) 

The thresholding strategy has to be adapted. 
Indeed, for the impainting task we need to have 
a very large initial threshold in order to have 
a very smooth image in the beginning and to 
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refine the details progressively. We chose an ex- 
ponentially decreasing threshold: 

A„ = A,nax(2'T';^^^^ - 1), n = 1, 2, • • • , iV,„ax, 

(24) 

where Amax = niax($ X). 



ficients of B at scale j and position k. We deter- 
mine the multi-resolution support by comparing 
\dj[k] — (i^-^^[fc]| with KCTj. 

We formulate the reconstruction problem as 
a convex constrained minimization problem: 



Algorithm 4 MS- VST + lUWT Denoising + 
Inpainting 

Require: data aq :— Y, mask 11, number of 
iterations iVmax, threshold k. 
Detection 



for j = 1 to J do 

Compute ttj and dj using 
Hard threshold |dj[fc]| with threshold naj 
and update A4. 
4: end for 
Estimation 



9: 
10 
11 



Initialize X'"' 0, Ao = A„ 



for n = to N„ 



1 do 



X = p+[x(") + *PM*^n(Y - x("))]. 

X("+l) ^ $ST[$7^X]. 



A„+l = Amax(2^ ^ 

end for 

Get the estimate A = X^^"'--"- 



Experiment 

We applied this method on simulated Fermi 
data where we masked the most luminous 
sources. 

The results are on Figure El The MS- VST 
+ lUWT + Inpainting method (Algorithm H]) 
interpolates the missing data very well. Indeed, 
the missing part can not be seen anymore in 
the inpainted map, which shows that the diffuse 
emission component has been correctly recon- 
structed. 



6. Source detection: denoising and 
background modeling 

6.1. Method 

In the case of Fermi data, the diffuse gamma- 
ray emission from the Milky Way, due to inter- 
action between cosmic rays and interstellar gas 
and radiation, makes a relatively intense back- 
ground. We have to extract this background in 
order to detect point sources. This diffuse in- 
terstellar emission can be modelled by a linear 
combination of gas templates and inverse comp- 
ton map. We can use such a background model 
and incorporate a background removal in our 
denoising algorithm. 

We note Y the data, B the background we 
want to remove, and df^ [k] the MS-VSTS coef- 



Argniin||*-'X||i,s.t. 
X > 0, 

V(j, fc) e X, (*^X),[fc] = (*^(Y - B)),[fc], 

(25) 



Then, the reconstruction algorithm scheme 
becomes: 

X = P+[X(") + *Pm*^(Y - B - X("))], (26) 
X(»+i) = $STa„ [*'^X]. (27) 

The algorithm is illustrated by the theo- 
retical study in Figure El We denoise Poisson 
data while separating a single source, which 
is a Gaussian of standard deviation equal to 
0.01, from a background, which is a sum of two 
Gaussians of standard deviation equal to 0.1 and 
0.01 respectively. 



Algorithm 5 MS-VSTS 
Background extraction 



lUWT Denoising 



Require: data cq := Y, background i?, 
number of iterations A^max, threshold n. 
Detection 



for j = 1 to J do 

Compute ttj and dj using ([7]). 

Hard threshold {dj[k] - rff'[fc]) with 



threshold and update M. 
4: end for 
Estimation 



6 
7 

8 
9 

10: 
11 



Initialize X'") = 0, Aq = 1. 
for n = to A'^inax — 1 do 

X = P+ [X(") + *Pm *^ (Y - B - X(") )] . 



= *STa„ [*'^X] 

A'm.x-("+l) 
Wmax-1 ■ 



An+l = 

end for 

Get the estimate A = X*^^""' 



Like Algorithm [1] Algorithm [H can be 
adapted to make multiresolution support adap- 
tation. 



6.2. Experiment 

We applied Algorithms [51 on simulated Fermi 
data. To test the efficiency of our method, we 
dete ct the sources with th e SExtractor rou- 
tine (jBertin fc Arnoutg|ll996( ). and compare the 
detected sources with the theoretical sources 
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catalog to get the number of true and false de- 
tections. Results are shown on Figures [7] and [HI 
The SExtractor method was applied on the first 
wavelet scale of the reconstructed map, with a 
detection threshold equal to 1. It has been cho- 
sen to optimise the number of true detections. 
SExtractor makes 593 true detections and 71 
false detections on the Fermi simulated map re- 
stored with Algorithm!!] among the 1000 sources 
of the simulation. On noisy data, many fluc- 
tuations due to Poisson noise are detected as 
sources by SExtractor, which leads to a big num- 
ber of false detections (more than 2000 in the 
case of Fermi data). 

6.2.1. Sensitivity to model errors 

As it is difficult to model the background 
precisely, it is important to study the sensi- 
tivity of the method to model errors. We add 
a stationary Gaussian noise to the background 
model, we compute the MS-VSTS + lUWT with 
threshold Scj on the simulated Fermi Poisson 
data with extraction of the noisy background, 
and we study the percent of true and false de- 
tections with respect to the total number of 
sources of the simulation and the signal-noise 
ratio (SNR(dS) = 201og( 

^signal I ^ noise 

)) versus 

the standard deviation of the Gaussian pertur- 
bation. Table [3] shows that, when the standard 
deviation of the noise on the background model 
becomes of the same range as the mean of the 
Poisson intensity distribution (Amean = 68.764), 
the number of false detections increases, the 
number of true detections decreases and the sig- 
nal noise ratio decreases. While the perturbation 
is not too strong (standard deviation < 10), the 
effect of the model error remains low. 



7. Conclusion 

This paper presented new methods for 
restoration of spherical data with noise follow- 
ing a Poisson distribution. A denoising method 
was proposed, which used a variance stabiliza- 
tion method and multiscale transforms on the 
sphere. Experiments have shown it is very ef- 
ficient for Fermi data denoising. Two spheri- 
cal multiscale transforms, the wavelet and the 
curvelets, were used. Then, we have proposed 
an extension of the denoising method in or- 
der to take into account missing data, and we 
have shown that this inpainting method could 
be a useful tool to estimate the diffuse emis- 
sion. Finally, we have introduced a new denois- 
ing method the sphere which takes into account 
a background model. The simulated data have 
shown that it is relatively robust to errors in the 
model, and can therefore be used for Fermi dif- 
fuse background modeling and source detection. 
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Table 1. Precomputed values of the variances aj of the wavelet coefficients. 



Wavelet scale j 


Value of Uj 


1 


0.484704 


2 


0.0552595 


3 


0.0236458 


4 


0.0114056 


5 


0.00567026 




Figure 2. Comparison of MS-VSTS with Anscombe + wavelet shrinkage on a single HEALPix 
face. Top Left ; Sources of varying intensity. Top Right : Sources of varying intensity with Poisson 
noise. Bottom Left : Poisson sources of varying intensity reconstructed with Anscombe + wavelet 
shrinkage. Bottom Right : Poisson sources of varying intensity reconstructed with MS-VSTS. 
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Table 2. Asymptotic values of the variances a-j^k of the curvelet coefficients. 



j 


I ^ 1 


I = 2 


I = 3 


I = 4 


1 


1.74550 


0.348175 






2 


0.230621 


0.248233 


0.196981 




3 


0.0548140 


0.0989918 


0.219056 




4 


0.0212912 


0.0417454 


0.0875663 


0.20375 


5 


0.00989616 


0.0158273 


0.0352021 


0.163248 



Simulated Fermi data without i 



Simulated Fermi data with. Poi 




Simulated Fermi data denoised with Anscombe + wavelet shrinkagi 



Simulated Fermi data denoised with IIS-VSTS +- Curvelets 





Simulated Fermi data denoised with MS-VSTS -+ lUWT Tvitli tliresliold 5sign: 



Simulated Fermi data denoised with MS-VSTS -+ lUWT ivith threshold Ssigm 




Figures. Top Left : Fermi simulated map without noise. Top Right : Fermi simulated map with 
Poisson noise. Middle Left : Fermi simulated map denoised with Anscombe VST + wavelet shrink- 
age. Middle Right : Fermi simulated map denoised with MS-VSTS + curvelets (Algorithm [3]). 
Bottom Left : Fermi simulated map denoised with MS-VSTS + lUWT (Algorithm [j) with thresh- 
old 5(Tj . Bottom Right : Fermi simulated map denoised with MS-VSTS + lUWT (Algorithm [IJ 
with threshold 3a j. Pictures are in logarithmic scale. 
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Figure 4. View of a single HEALPix face from the results of Figure [3l Top Left : Fermi simulated 
map without noise. Top Right: Fermi simulated map with Poisson noise. Middle Left: Fermi simu- 
lated map denoised with Anscombe VST + wavelet shrinkage. Middle Right : Fermi simulated map 
denoised with MS- VSTS + curvelets (Algorithm [3)) . Bottom Left : Fermi simulated map denoised 
with MS- VSTS + lUWT (Algorithm |T]) with threshold 5cr,. Bottom Right : Fermi simulated map 
alanoised with MS- VSTS + lUWT (Algorithm H]) with threshold 3o-j . Pictures are in logarithmic 
scale. 
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Simulated Fermi map with missing data Simulated Fermi data deaoised and iapainted by MS— VSTS+HTWT 




^ 1^ 



Figure 5. MS-VSTS - Inpainting. Left : Fermi simulated map with Poisson noise and the most 
luminous sources masked. Right: Fermi simulated map denoised and inpainted with wavelets 
(Algorithm H]) . Pictures are in logarithmic scale. 



Table 3. Percent of true and false detection and signal- noise ratio versus the standard deviation 
of the Gaussian noise on the background model. 



Model error std dev 


% of true detect 


% of false detect 


SNR (dB) 





59.3% 


7.1% 


23.8 


10 


57.0% 


11.0% 


23.2 


20 


53.2% 


18.9% 


22.6 


30 


49.1% 


43.5% 


21.7 


40 


42.3% 


44.3% 


21.0 


50 


34.9% 


39.0% 


20.3 


60 


30.3% 


37.5% 


19.5 


70 


25.0% 


34.6% 


18.9 


80 


23.0% 


28.5% 


18.7 


90 


23.6% 


27.1% 


18.3 
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Figure 6. Theoretical testing for MS-VSTS + lUWT denoising + background removal algorithm 
(Algorithm El . View on a single HEALPix face. Top Left : Simulated background : sum of two 
Gaussians of standard deviation equal to 0.1 and 0.01 respectively. Top Right : Simulated source: 
Gaussian of standard deviation equal to 0.01. Bottom Left : Simulated poisson data. Bottom Right : 
Image denoised with MS-VSTS + lUWT and background removal. 
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Milky-Way background model Simulated Femi point sources 




Figure 7. Top Left : Simulated background model. Top Right : Simulated Gamma Ray sources. 
Middle Left : Simulated Fermi data with Poisson noise. Middle Right: Reconstructed Gamma Ray 
Sources wfth MS-VSTS + lUWT + background removal (Algorithm^ with threshold Scr^ . Bottom : 
Reconstructed Gamma Ray Sources with MS-VSTS + lUWT + background removal (Algorithm [5]) 
with threshold 3a j. Pictures are in logarithmic scale. 
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Figures. View of a single HEALPix face from the results of Figure [71 Top Left : Simulated back- 
ground model. Top Right : Simulated Gamma Ray sources. Middle Left : Simulated Fermi data with 
Poisson noise. Middle Right: Reconstructed Gamma Ray Sources with MS-VSTS + lUWT + back- 
ground removal (Algorithm [S]) with threshold 5a j. Bottom: Reconstructed Gamma Ray Sources 
with MS-VSTS + lUWT + background removal (Algorithm O with threshold 3a j. Pictures are in 
Icggarithmic scale. 



